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Abstract 

In this paper we present continuous age- and space-structured models and numerical 
computations of Proteus mirabilis swarm-colony development. We base the mathemat- 
ical representation of the cell-cycle dynamics of Proteus mirabilis on those developed 
by Esipov and Shapiro, which are the best understood aspects of the system, and we 
make minimum assumptions about less-understood mechanisms, such as precise forms 
of the spatial diffusion. The models in this paper have explicit age-structure and, 
when solved numerically, display both the temporal and spatial regularity seen in ex- 
periments, whereas the Esipov and Shapiro model, when solved accurately, shows only 
the temporal regularity. 

The composite hyperbolic-parabolic partial differential equations used to model 
Proteus mirabilis swarm-colony development are relevant to other biological systems 
where the spatial dynamics depend on local physiological structure. We use compu- 
tational methods designed for such systems, with known convergence properties, to 
obtain the numerical results presented in this paper. 

Key words: Proteus mirabilis, sweirm colony, age structure, space structure, par- 
tial differential equations . 

1 Introduction 

In this paper we present continuous age- and space-structured models and numerical compu- 
tations of Proteus mirabilis swarm-colony development. Proteus mirabilis has received treat- 
ment in the mathematical biology literature due to the strikingly regular spatial-temporal 
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patterns swarm colonies make on an agar surface ^1 HI] . Not only do the swarm colonies 
form regularly spaced concentric terraces in a bulls-eye pattern, but the total time taken to 
form a ring is invariant under changes in the agar or glucose concentration in the substrate. 

The invariance of the total cycle time rules out chemotaxis as an explanation for the 
ring formation, as do biological observations of Proteus responses to chemical stimuli |49j . 
Esipov and Shapiro ^Hl saw that the regularity of the ring-formation cycle is a result of the 
regularity of the differentiation-dedifferentiation cycle of individual bacteria. The resulting 
model was an equation that is first-order hyperbolic in the age variable (also called a kinetic 
or transport term) and nonlinear parabolic in space. 

This mathematical explanation for swarm-colony behavior necessarily differs from chemo- 
taxis models, such as those for swarm-colony formation in E. coli, due to the temporal 
invariance of the Proteus system j4Uj . 

Proteus has been suggested as a model system for problems in morphogenesis and other 
complex self-organizing systems (THl EHl Hll , but it is also an example of a biological system 
where the local kinetics - the age structure - determine the spatial dynamics of the system. 
Dictostelium has received much treatment as a model system for other problems in mor- 
phogenesis |22|- The study of Proteus can have similar relevance to systems such as tumor 
growth 13 [T31 EHl and natural forests [3 |2Zl 123 EOl , where physiological structure, such as 
age, size or activation levels in tumor and immune cells, can play an important role in the 
spatial behavior of the system. Bertuzzi and Gandolfi presented a cancer tumor-cord 
model with both space and size structure with some similarity to the kinetics in the Proteus 
models, although they did not present results for their full model. 

We base the models of Proteus swarm-colony development in this paper on models pre- 
sented in j^. Due to the phonetic and conceptual similarities of the terms "swimmer" and 
"swarmer", we depart from standard terminology and use the term "dividing cell" instead 
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of "swimmer cell" . This change in terminology emphasizes both the nonmotile character of 
the dividing cells on an agar surface and the fact that the swarmer cells are filaments^. 

The models in jl6j indicated that experimental study of sharp dedifferentiation in Pro- 
teus might yield insight into their periodic swarm-colony development. Dedifferentiation, also 
called septation, is the process in which multinuclear swarmer cells subdivide into mononu- 
clear dividing cells. This experimental study remains important for the models in this paper, 
but perhaps more critical are experiments that determine how the lag phase in differentia- 
tion - the process in which mitosis of dividing cells results in a swarmer cell rather than two 
dividing cells - affects swarm-colony development, and ideally give the functional forms of 
this lag. This lag in the onset of differentiation was identified in ^^-i and is an important 
part of the models in this paper and of those in 39 , but is not used in |T^ . 

We begin the paper with a description of Proteus swarm-colony development and aspects 
of the relevant biology, and then discuss differences between the models in this paper and 
prior work. When referring to the Esipov-Shapiro system, we mean a modified form of the 
system formally presented in JH] • Without these changes, the Esipov-Shapiro system has no 
swarming behavior at all and is thus not really a model of the swarm-colony development. 
In Appendix |XJ we elaborate upon the changes and why they might reflect the intended 
mechanisms of Esipov and Shapiro. Again, without these changes, we can have no useful 
discussion of their model. 

The most significant difference between the models in this paper and those in JB] is that 
the models in this paper reproduce the spatial regularity seen in Proteus swarm colony devel- 
opment, whereas, when solved accurately, the Esipov-Shapiro system, even when modified 

with the maximum benefit of doubt, does not. The damping in terrace width is discussed in 

^ James Shapiro at the University of Chicago suggested this terminology when discussing the behavior of 
Proteus on an agar surface. 
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more detail in Appendix ^ 

In subsequent sections, we present our mathematical models, the numerical methodology 
used to solve the age-, space- and time-dependent systems, and the results from the numerical 
computations. 

2 Proteus mirabilis Swarm- Colony Development 

When inoculated onto an agar surface, Proteus begin formation of strikingly regular spatio- 
temporal patterns that begin with three initial phases: a lag phase, a first swarming phase, 
and a first consolidation phase; these form what appear to be the center circle and first 
ring of a bulls-eye. This is followed by repeating cycles of ring, or terrace, formation that 
consist of a swarming phase followed by a consolidation phase. These subsequent terraces 
are equal in width and time of formation. Perhaps most interestingly, the time of terrace 
formation is invariant under changes in the glucose or agar concentration of the substrate, 
although changes in agar concentration do affect the ratio of time spent swarming versus 
consolidating, and changes in glucose increase terrace width. Terrace formation does vary 
by temperature jSj. 

Broth cultures of Proteus consist of mononuclear cells (often called "swimmer" cells, but 
in this paper called "dividing" cells) about 0.6 fim wide and 1-2 fim long with short flagella. 
On an agar surface these cells lengthen somewhat to about 0.8 fim wide and 2-4 /xm long. 
After a period of time^, some cells near the perimeter begin to lengthen into multinuclear 
filament cells (called "swarmer cells") that reach sizes of 0.7 fim wide and 20-80 /im long. The 
swarmer cells have longer fiagella than the dividing cells. Moreover, dividing cells normally 

have between 1 to 10 fiagella, whereas swarmer cells have between 500 to 5000, a 10-fold 

^This time-lag before swarmer cell production begins is incorporated into the models in this paper and 
in j39| as a lag in density, but is not used in 16 . 



4 



increase in the number of flagella per cubic micrometer of bacterial volume jlH] • The process 
in which dividing cells become swarmer cells is called "differentiation." Differentiation only 
occurs above one dividing-cells density and below another. 

It is the swarmer cells that move across the agar during the swarming phase. The 
movement of swarmer cells on the agar surface is a complex process. The swarmers use their 
elongated flagella to move in fluid extracted by them from the agar gel. They do so in groups, 
or "rafts" ^H]- Consequently, the unit of random motion in any diffusion approximation of 
Proteus movement is a large and complicated entity. 

When the multinuclear swarmer cells approach a maximum size, they rapidly break down 
into single- nucleus dividing cells. This process is called "dedifferentiation." Since mitosis 
results in exponential growth, the size of a swarmer cell is an exponential weight times its 
age. Thus, dedifferentiation occurs when swarmers reach a maximum age. 

Varying glucose concentration greatly alters biomass production without greatly affecting 
the spatial or temporal aspects of the swarm and the consolidation cycle of terrace formation. 
Thus, the cycle is not affected by biomass. Changes in agar concentration do not change 
the total time of terrace formation, but higher agar concentrations shorten the time spent 
swarming and the width of the terraces, and lengthen the time spent in consolidation |TH]. 

3 Differences and Commonalities with Previous Mod- 
els 

Models presented in this paper differ from models presented in |I6| in several respects. We 
use a radially symmetric geometry with no-fiux boundary conditions, as would be the case 
on a petri dish. Models used in (THj use a linear geometry. We use a diffusion approximation 
which arises from isotropic random motion instead of Fickian diffusion ^1 E]- An 
additional difference between the models in this paper and in is that we present results 
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for nondimensionalized models. 

The most substantial differences between the models presented here and in ^H] are in the 
different forms of the diffusivity, D, and the differentiation rate to swarmer cells, ^ (as will be 
discussed in Sectional our numerical treatments also differ, with nontrivial consequences.) 
The diffusivity in relied on a memory field: the diffusivity changes depending if the cells 
are in a collective state of motion or at rest. This is not clearly supported biologically and, 
invoking Occam's Razor, is not used in the models in this paper, since it is not needed. Also, 
the diffusivity in (TE] depends on the dividing-cell population. 

The diffusivity in this paper depends only on swarmer biomass. Since we use swarmer 
biomass rather than swarmer numerical density in the diffusivity, the exponential weighting 
from mitosis that describes size means that younger swarmer cells do not contribute as much 
as older swarmer cells. This weighting also models the fact that each nucleus of a swarmer 
cell dedifferentiates into a dividing cell. Because of this weighting, age-structured models 
of Proteus swarm-colony development presented in this paper can also be thought of as 
size-structured. 

The second and most important difference, the form of ^, is based on incorporating the 
property that dividing cells do not form swarmer cells until a minimum population density 
is reached, a mechanism discussed in jlH] and introduced in models in [;39j, but not used in 
jl6j . This threshold is represented mathematically by setting ^ to zero for low dividing-cell 
densities. 

We focus on three features of Proteus swarm-colony development: cycle regularity, ring- 
width uniformity, and control of the ratio of swarm to consolidation time. A regular cycle 
does occur in the Esipov-Shapiro system, and these authors are able to control the ratio of 
swarm time to consolidation time. However, the terrace width decreases from cycle to cycle, 
so the front velocity does not remain periodic, as shown in Appendix 1X1 
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In jin] it was concluded that models that used a sharp age of septation from swarmer 
filament to dividing cells gave periodic front dynamics, whereas a constant dedifferentiation 
rate did not. 

Medvedev, Kaper, and Kopell [Sn] developed a reaction-diffusion model of Proteus swarm- 
colony development based on averaging over the age variable. The model relied on a specific 
nonlinear form of the diffusivity to produce periodic front dynamics. In |39j , a constant ded- 
ifferentiation rate was used rather than a sharp age of septation. To obtain a periodic front 
velocity, the models in also introduced a piecewise constant differentiation rate with a 
lag phase - a biological mechanism mentioned in j49j and used in the models in this paper - 
and a diffusivity dependent on the ratio of swarmer-cell to dividing-cell biomass. The use of 
a constant dedifferentiation rate allowed the age-structured diffusion system to be reduced 
to a reaction-diffusion system. This reduction has advantages both numerically and analyt- 
ically. Moreover, a reaction-diffusion equation with periodic front dynamics is an interesting 
mathematical object in its own right. However, removing explicit age eliminated a mecha- 
nism for controlling the ratio of time spent swarming to time spent in consolidation without 
changing the total cycle time. No mechanism within the reaction-diffusion framework was 
suggested to replace this. 

The biological mechanisms that are best understood are emphasized in the models in 
this paper. In particular, we maintain the explicit age structure. Due to experimental 
results'^ which postdate previous work I39j. we are more confident in assuming a sharp 
age of septation. New numerical methods for handling age-structured population models 
with nonlinear diffusion have been developed so that computations that incorporate the 

known properties of Proteus swarm-colony development can be run efficiently |3 El • 

•^Preliminary experiments conducted by James Shapiro's laboratory at the University of Chicago suggest 
that Proteus mirabilis that have been genetically modified to no longer have a sharp septation age lose 
periodicity in their swarm-colony development. 
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Because less is known about the specific forms of the diffusion and diffusivity, we attempt 
to use the simplest form of diffusivity suitable for this system. Rather than using a relatively 
complicated fractional form involving both swarmer and dividing cells, we use a diffusivity 
that is proportional to swarmer biomass above a threshold, reflecting the need for swarmers 
to group together in order to move. This is the simplest form that is biologically justified 
and that results in both swarming and consolidation. This density dependence on swarmers, 
due to "raft" building, is not unique to the models in this paper and is also used in jTH] and 
|39j . The diffusivity we use, since ^ is zero above some threshold, also results in an absence 
of swarmers in the interior of the colony - the primary justification for the dependence of 
the diffusivity on dividing cells in the other models. 

The different model mechanisms and results for those in this paper, Esipov-Shapiro, 
and Medvedev et al. are summarized in Table ^ The models in this paper are unique in 
reproducing three major aspects of Proteus swarm-colony development: temporal regularity 
of overall cycle time, spatial regularity, and control of the ratio of swarm time to consolidation 
time. Moreover, by retaining the explicit age structure in the mathematics, the models 
reproduce the colony behavior by relying on only two main mechanisms. 

4 The Model 

The models in this paper and in ^B] are based on a mathematical framework with some 
history behind it. Skellam 46^ considered the effects of diffusion on populations in his classic 
work of 1951. Sharpe and Lotka in 1911 and McKendrick in 1926 considered population 
models with linear age structure [48J. More recently, Gurtin and MacCamy [201 considered 
models with nonlinear age structure. Rotenberg [ISj and Gurtin [THj posed models dependent 
on both age and space. Gurtin and MacCamy [21] differentiated between two kinds of 
diffusion in these models: diffusion due to random dispersal, and diffusion toward an area 
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Mechanism \ Model 


Ayati 


Esipov-Shapiro 


Medvedev et al. 


Dependence of diffusivity 








on swarmer density 








Dependence of diffusivity 








on dividing-cell density 








Memory field 








Lag phase in ^ 








Explicit age structure 




v/ 




Behavior \ Model 


Ayati 


Esipov-Shapiro 


Medvedev et al. 


Regularity of swarm and 








consohdation cycle 


a/ 






Regularly spaced 








concentric terraces 








Control of ratio of swarm 








to consolidation time 









Table 1: Mechanisms and behaviors for different models of Proteus swarm-colony develop- 
ment. The models in this paper are unique in reproducing three major aspects of the colony 
development, and, by retaining the explicit age structure, reproduce the colony behavior 
with only two main mechanisms. 
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of less crowding. Existence and uniqueness results can be found for various forms of these 
models in Busenberg and lannelli [Hj, di Blasio di Blasio and Lamberti J3], Langlais 
|33j . and MacCamy Further analysis has been done by several authors 1^ ISH I38j . 

For our model, we assume the colony is radially symmetric. The variables r, a, and t 
represent radius in two-dimensional space, age, and time, respectively. The function u{r, a, t) 
represents the swarmer-cell population density at radius r, age a, and time t. The functions 
p{r,t) and v{r,t) represent the biomass density of swarmer cells and the dividing-cell popu- 
lation density, respectively, at radius r and time t. The dividing-cell population density, v, is 
measured in population-density units (pdu) with pdu = Gc/cm^, where Gc denotes one bil- 
lion cells. The swarmer-cell population density, u, is measured in pdu/hr. The swarmer-cell 
biomass density, p, is measured in mg/cm^. 

The dimensional model consists of the age-structured nonlinear diffusion system 

dtu + daU = -dr{rdr{D{p)u)) — fi{a)u, < r < r^, a > 0, t > 0, (4.1a) 



. l-e(t^) 

OtV = V 

T 



poo 

+ fx{a)ue''/^ da,0 <r <rm,t > 0, (4.1b) 
Jo 

u{r, 0, t) = ^^v{r, t), < r < r^, t>Q (4.1c) 

T 



with boundary and initial conditions 

dr{D{p{r^,t))u{rm,a,t)) =0, a>0,t>0, (4.1d) 

u{r, a, 0) = 0, <r <rm,a>0, (4.1e) 

f (r, 0) = Vo{r), <r <rm. (4. If) 

The parameter r is the time it takes a cell to subdivide. A typical value of r is 1.5 hours. 
The total motile swarmer cell biomass is given by 

POO 

p(r,t)=mo u{r,a,t)e"'''' da, 0<r<rm,t>0, (4.1g) 
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where mo is measured in mg/Gc. A dividing cell has a volume of approximately 2 x 10~^^m^ 
so that mo ~ 2 for Proteus mirabilis. 

Proteus move through a process of raft building that requires two things: sufficient 
maturity in swarmer cells to contribute to raft building, and a sufficient biomass of mature 
cells to form the rafts. The lower limit of integration, a^^^, is the minimum age when a 
swarmer cell is sufficiently large, with sufficiently long flagella, to contribute to motion on 
the agar surface. The parameter a^^^ is comparable to the parameter 6^^^, the lower limit of 
integration in the total swarmer numerical density, in ^Hl- Each parameter controls the ratio 
of swarm time to consolidation time without changing the total time of ring formation. The 
parameter a^^i^ is related to agar concentration. The higher the concentration, the drier the 
surface, raising the value of a„,i„. Higher agar concentration, and thus higher a^i^, shortens 
the swarming phase and lengthens the consolidation phase in terrace formation. 

The diffusivity has the form 



The parameter p^^^ is the minimum biomass needed for swarmers to build rafts capable of 
moving on the agar surface, and is roughly equivalent to P^ j^^^j^ in jTHj. 

We use a differentiation function with a lag phase that is a piecewise cubic with 
support of length 2vw- 



The interval [vc — Vw,Vc + v^] is the swarmer-cell production window. The parameter Vc 
represents the lag phase in swarmer-cell production We use a compactly supported 

cubic function because it is smooth. The important aspects of the functional form of ^ are 
that it is zero above and below certain thresholds. Appendix iBl shows that the nature of the 



D{r, t) = Dq max 



{{p{r,t) -Pn,in),0}. 



(4.1h) 




(4.1i) 
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results do not depend strongly on the shape of the curve. The shape of the curve does help 
with numerical issues concerning the degenerate diffusion. The parameter Vc has no analog 
in Unj. 

The dedifferentiation modulus, with a spread parameter a, is 

/i(a) = — exp ((a - a^,jV (2cr^)) + —H{a), (4.1j) 
0" cr 

where 

^ ^ [0, a < 

is the Heaviside function. This dedifferentiation modulus represents a situation where the 
probability of dedifferentiation is low for young swarmers, increases as they approach a^^^ in 
age, and remains high afterwords. Experiments indicate that this type of dedifferentiation 
is dominant for Proteus mirabilis, and that the transition, from virtually no probability of 
dedifferentiation to an extremely high probability of dedifferentiation as a function of age, 
occurs very rapidly near the critical age, a^^^. Preliminary results from James Shapiro's 
laboratory at the University of Chicago show that different strains of Proteus with less sharp 
dedifferentiation have oscillatory, but not strictly periodic, front velocities. 
We take the initial condition to be 

Vo{r) = { -l^V-; -V-)' + 0' ° (4.1k) 




0, r > tq. 

Again, it is mass, not shape, that matters qualitatively. Shape effects numerical efficiency. 

In the limit a ^ 0, which represents a sharp age at which all dedifferentiation occurs, we 

have a system on a bounded age domain 

dtU + daU = ^dr{rdr{D{p)u)), < r < rm,0 < a < a^,^,t > 0, (4.2a) 

dtV = ^—^v + u{r, a^,., t)e"— < r < r„, t > 0, (4.2b) 
r 

p{r,t) = mo ue"/^ da, <r <r^, t>0 (4.2c) 



a„ 
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with conditions 

u{r,0,t) = ^v{r,t), <r <r^, t>0, (4.2d) 
r 

dr{D{p{rra,t))u{rm,a,t)) = 0, < a < a^,,, t > 0, (4.2e) 

u{r, a, 0) = 0, < r < r^, < a < a_,, (4.2f) 

i;(r,0) = wol?"), < r < r^, (4.2g) 

4.1 A Nondimensional Model 

We nondimensionalize Equations ()4.1a|l - ()4.1kj) . We scale age and time by r and space by 
and introduce the independent variables 

r = r/vm, a = a/r, t = t/r. 
We scale the dependent variables, 

U{r,a,i) = u{r,a,t)/{v^/r), V{f,i) = v{r,t)/v^, P = p/{moV^). 
We use the dimensionless parameters 

^max (^inax/'T'l ^min ^ <J j T . 

We note that in this nondimensionalization, the length of the swarmer-cell production 
window (2f^) and swarmer-cell swarm threshold (pmin) are related. It seems intuitive that 
halving one is equivalent to doubling the other. The effect of a lag in swarmer cell production 
(fc) is isolated from the effects of Vw and p^^^. 



13 



For (T 7^ 0, we obtain the system 

diU + daU = ^dr{fdr{DiP)U)) - il{d)U, < r < 1, a > 0, i > 0, (4.3a) 

/•oo 

dfV = {l-i{V))V+ il{a)Ue^ da, <r <1, i> 0, (4.3b) 

Jo 

P{f, i) = / U{f, d, f)e" dd, < r < 1, £ > (4.3c) 

Omin 



with conditions 



C/(f, 0, i) = |(F)F(f, £), < f < 1, £ > 0, (4.3d) 

dr{D{P{l,t))U{l,d,i)) ^0, a>0,£>0, (4.3e) 

C/(f, a, 0) = 0, < r < 1, a > (4.3f) 

l^(f, 0) = VQ{r^r)/v^, 0<r <rm, (4.3g) 



where 



D{f,i) = Domax{(P(f,t)-^.nin),0}, (4.3h) 

= Uo(2|i^-KP-3(i--K)^ + i), |v-v;|<i, 

^ ' 1^ 0, otherwise, ^ ' ' 

Ma) = ^ exp((a-a_)V(2<7^)) + ^i/(a). (4.3j) 

In the hmit it — > we have the system 

dfU + daU = ^dr{fdf{D{P)U)), 0<f<l, 0<a< a^,., £ > 0, (4.4a) 

diV = (1 - e(y))V^ + Uir, a_, t>"— , < f < 1, t > 0, (4.4b) 

P(f, £) = / C/e" da, < f < 1, £ > 0, (4.4c) 
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with conditions 



U{f,Q,i)=i{V)V{f,i) 



< f < 1, t > 0, 



(4.4d) 



9,(D(P(l,t))f/(l,a,t)) =0 



< d < a, 



'max 5 



t > 0, 



(4.4e) 



f/(f,a,0) = 0, 



0<f<l, 0<a<d, 



'max 7 



(4.4f) 



V{f,0) = Vo{rmf)/v, 



< f < 1. 



(4.4g) 



In this paper we will examine only the limiting case of a — > 0, which corresponds to the 
case of a sharp age of septation, "Model A", in [16j. The models used in the computations 
in this paper are nondimensional, whereas dimensional models were used in |16j . 

As mentioned earlier, there are two major differences between the models in this paper 
and those in J^l- First, the diffusivity is of much simpler form, with no memory field, and 
hence no parameter for the upper threshold for diffusion in a memory field (P^^„). Second, 
there is a lag phase in the dedifferentiation function, ^. Also, we use biomass, not numerical 
density (which was how the model was stated in [TBlj) as the relevant measure of the totality 
of swarmers. Other differences include radial symmetry, and diffusion based on isotropic 
random motion rather than Fickian diffusion. 

When solved accurately, these differences yield the temporal regularity seen in Proteus 
experiments, as well as the regularly spaced concentric terraces, something not seen in ac- 
curate computations of the Esipov-Shapiro system. 

5 Numerical Methodology 

There has been much investigation into numerical methods for solving models with just age 



structure HZl 1^21 IHSl HZ]- Kim Kim and Park ^26J, and Lopez and Trigiante 



developed methods for age-structured populations that undergo random dispersal in space. 
All these methods involve uniform time and age discretizations, with the age step chosen 
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to equal the time step. These methods discretized along characteristics, but they did so 
simultaneously in age and time and thus imposed the often crippling constraint that the 
time and age steps must be both constant and equal. The difficulty with this approach is 
twofold. First, the use of constant age and time steps prevents adaptivity of the discretization 
in age or, especially, time. Second, and more important, the coupling of the age and time 
meshes can cause great losses of efficiency, since only rarely will the dynamics in time be on 
the same scale as the dynamics in age. This is particularly the case when space is involved, 
since sharp moving fronts can require small time steps, whereas the behavior in the age 
variable can remain relatively smooth. 

The methods used in this paper |21 El use a moving age discretization that allows for 
nonuniform age and time discretizations. The age step need not equal the time step. Instead, 
the positions of the age nodes are adjusted by the time step. The methods preserve the 
important fact that age and time advance together. The age discretizations presented in 
123 123 ISHI can be viewed as special cases of the methods presented in [21 El by setting the 
time and age meshes to be constant and equal and using a backward Euler discretization in 
time and a piecewise constant finite element space in age. 

The importance of allowing different age and time discretizations is somewhat illustrated 
by the application of the methods of de Roos ^U]. These methods have found use in the 
study of ecological systems such as Daphnia (see [TT] and the references therein) as well as 
in theoretical population biology (23 112] ■ The methods of de Roos are formulated for the 
case of time and a variable representing some sort of physiological structure, most simply 
age, and involve moving the age nodes along characteristics. They have not been formulated 
for explicit space. The representation of the approximate solution is probabilistic and not 
functional, and birth and death are handled differently than in the methods in [21 E] . Even 
so, it would be interesting to know if an energy analysis could provide a framework for the 
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convergence analysis sought in [T2]. The main effect of de Roos's methods is to separate 
the age and time discretizations, thereby yielding an approximation that is dispersion free 
in age, in order to provide a method that works in practice. 

We emphasize the importance of reliable numerics. In the systems presented in this 
paper and in ^^I; the regularity in terrace formation is due to the regularity in the aging 
of the swarmer cells. A coarse, stage-structured numerical approximation to the continuous 
aging term, such as that used in fW], can constitute a qualitatively different swarmer-cell 
life cycle than that of the original continuous model. Also, the degenerate diffusion can be 
altered by the choice of the time step. The numerical computation may then show periodic 
front dynamics that are not obtained by an accurate solution of the original continuous 
equations, but are induced by a regularity in the numerics. Appendix 1X1 details how accurate 
computations change the behavior of the Esipov- Shapiro system from what is presented in 

m- 

Proteus mirahilis swarm-colony development is important as a model system for other 
problems in morphology, complex self-organizing systems, and where spatial behavior is a 
manifestation of local kinetics. Not only are the choices made in the mathematical models 
important in this respect, but so is the computational framework used to solve them. One 
goal of this paper is to illustrate the necessity of robust and accurate new numerical methods 
in the modeling and solution of such systems. 

The results of the numerical computations presented in Section IHl were obtained for the 
various parameter sets by using the age and space discretizations described in [21 E] with 
time integration done using step-doubling extrapolation. Step-doubling is a time-stepping 
algorithm, going back to at least Gear for ordinary differential equations, based on taking 
one step over a time interval to obtain one approximate solution for that interval, and then 
taking two half-steps over the same interval to get a second approximate solution. The two 
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solutions can then be compared for adaptive step-size control. Moreover, the two solutions 
can be combined to obtain a likely second-order approximation in time j3]. The use of a 
second-order adaptive method with error checking for these problems is important, since the 
front dynamics of a degenerate system can depend on the time step. 

A uniform spatial discretization of size Ax = 1/300, and a uniform age discretization of 
Aa = 1/40, using piecewise constant basis functions, was sufficient to solve the system within 
a relative error in the L^-norm of less than 1%. A uniform age discretization in the context 
of a moving grid method means that all but the first and last age intervals are constant in 
length, and that a new age interval is introduced at the birth boundary when the old birth 
interval reaches Aa in length. Convergence in time was obtained by adjusting a tolerance 
parameter for the adaptive time-stepping in the step-doubling algorithm so that the relative 
error in the L^-norm was also less than 1%. 

6 Results 

We examine the response of the model to changes in the various parameters in the following 
tables. The spatial domain of the problems is f G [0, 1], the temporal domain is truncated 
to the finite domain t G [0,35], and the age domain is d G [0,dn,ax], which varies depending 
on the age of sharp septation in swarmers. 

The quantity T is the total time it takes for a swarm ring to form, S is the time spent in 
a swarm phase, C is the time spent in a consolidation phase, and R is the width of a swarm 
ring. Except for the parameter under study, the other parameters are set to a^^^ = 2.67, 
Vc = 8, Dq = 2 X 10^'^, P^i^ = 0.5, ^0 = 0.5, and a^i^ = 0. In each of the tables, the 
default value of the parameter will be shown in bold type. Figure Q shows the results of the 
model for the base parameters. The top graph is a three dimensional plot of the physical 
colony viewed from above and the bottom graph is a plot of colony radius versus time that 
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Figure 1: Computed Proteus mirabilis swarm colony with regularly spaced concentric terraces 
and periodic front velocity. The parameters for this swarm colony are d^^^ = 2.67, Vc = 8, 
Dq — 2 X 10~^, Pniin = 0.5, ^0 = 0.5, and a^^^ — 0. The top figure is a 3D plot of the 
swarmer- and dividing-cell biomass viewed from directly above with a combination of ambient 
light, diffuse reflection, and specular reflection. The bottom figure shows the radius of the 
swarm colony as a function of time. Areas where the function has zero slope correspond to 
consolidation periods in the swarm-colony development. In dimensionless units, total time 
for the formation of one ring is T = 4.1, the time spent swarming is S — 2.2, the time spent 
in consolidation is C = 1.9, and the width of one ring is i? = 0.153. 
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S/C 


R 


R/S 


0.0 


4.1 


2.2 


1.9 


1.16 


0.153 


0.0695 


0.5 


4.1 


2.0 


2.1 


0.95 


0.150 


0.0750 


1.0 


4.1 


1.7 


2.4 


0.71 


0.140 


0.0824 


2.0 


4.1 


1.3 


2.8 


0.46 


0.090 


0.0692 



Table 2: Ring formation function of a, 





T 


S 


C 


S/C 


R 


R/S 


0.001 


4.05 


2.15 


1.9 


1.13 


0.107 


0.0498 


0.002 


4.1 


2.2 


1.9 


1.16 


0.153 


0.0695 


0.003 


4.1 


2.3 


1.8 


1.28 


0.183 


0.0796 


0.004 


4.1 


2.3 


1.8 


1.28 


0.214 


0.0930 



Table 3: Ring formation function of Dq. 



illustrates the periodic front dynamics of Proteus swarm-colony development. 

Table El exhibits ring formation as a function of dmin, the minimum nondimensional age for 
a swarmer cell to contribute to swarming. The parameter a^i^ depends on agar concentration. 
As agar concentration increases, so does the difficulty of moving on the substrate, and hence 
dniin increases as well. As increases, 5* and R decrease, C increases, and T stays the 
same. This is the response to increased agar concentrations seen in experiments; as agar 
concentration increases, swarm time and terrace width decrease, consolidation time increases, 
whereas total cycle time stays the same. 

This invariance of the total time for ring formation, T, is further illustrated in Figure El 
A parameter similar to dmin controls the ratio of swarm time to consolidation time in |16j . 
Having averaged out the age variable, no such control is shown in [HH] . 

Table 01 exhibits ring formation as a function of Dq, the constant of diffusion. The 
parameter Dq reflects changes in glucose concentration. The temporal metrics, T, 5, and 
C, change little, whereas there is signiflcant change in the width of the terraces, R. This 
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Figure 2: Swarm-colony front dynamics as a function of d^i^. Shown are curves for the cases 
of a^in = 0.0, d,„in = 1.0, and a„,i„ = 2.0. The other parameters for these swarm colonies are 
d^^^ = 2.67, Vc = 8, Dq = 2 X 10~^, P„^i„ = 0.5, and = 0.5. The curves have been reahgned 
in time for better comparison of individual ring formation. Although total cycle time does 
not change, the ratio of time spent swarming to time spent in consolidation decreases as the 
minimum nondimensionalized age to contribute to swarming, a^i„, increases. 

corresponds to the experimental results illustrated in Figure 14 in which shows minor 
variation in the total cycle time and the ratio of swarm to consolidation time, while showing 
a significant increase in terrace width. 

Table lU exhibits ring formation as a function of Vc, the nondimensionalized lag in the 
swarmer cell production, ^. As Vc increases, we see an expected increase in the consolidation 
time, C. The time spent swarming, S, is nearly invariant. This leads to an increase in the 
overall cycle time, T, as Vc increases. The increased time building up swarmer rafts results 
in an increase in the front velocity, R/S, and thus in terrace width, R, as well. 

The presence of a lag in swarmer-cell production was observed in jiH] • However, the length 
of the lag, much less its explicit density dependence, is an open area of experimentation. The 
biological and chemical factors that underlie Vc are also unknown. 
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2.1 
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1.75 
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4 


3.6 


2.2 


1.4 


1.57 
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0.0605 


8 


4.1 


2.2 
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1.16 


0.153 


0.0695 


16 


4.5 
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2.3 


0.96 


0.157 


0.0714 



Table 4: Ring formation as a function of Vc- 





T 


S 


C 


5/C 




R/S 


0.25 


4.7 


1.65 


3.05 


0.54 


0.067 


0.0406 


0.50 


4.1 


2.2 


1.9 


1.16 


0.153 


0.0695 


0.75 


3.7 


2.6 


1.1 


2.40 


0.230 


0.0885 


1.00 


N/A 


N/A 


N/A 


N/A 


N/A 


N/A 



Table 5: Ring formation as a function of ^o- 



Table El exhibits ring formation as a function of the constant of the differentiation 
ratio, ^. The parameter is clearly greater than zero, otherwise dividing cells would never 
beget swarmers, and is thought to be less than one, as suggested in jlH], "some of the cells 
near the perimeter of the colony begin to undergo a dramatic morphological change to highly 
elongated forms. . . " Results for = 1 show the loss of a swarming consolidation cycle and 
a nearly constant rate of colony expansion. 

As .^0 increases, the consequent increase in swarmer-cell production results in longer 
time spent swarming, S", and terrace width, _R, but a decrease in total cycle time, T, and 
consolidation time, C. 

Table El exhibits ring formation as a function of P^in, the minimum swarmer biomass 
for swarming to begin. The parameter P^^^ reflects the need for swarmers to build "rafts" 
with other swarmers in order to move. "Movement of individual cells seems to be retarded 
whereas cells moving in larger groups do so much more effectively" [IHj. Rafts are able to 
form in sufficient numbers only when there is a minimum swarmer biomass present. The 
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N/A 


N/A 


0.186 


N/A 
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4.2 
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1.4 


2.00 


0.170 


0.0607 


0.3 


4.15 


2.6 


1.55 


1.68 


0.163 


0.0627 


0.5 


4.1 


2.2 


1.9 


1.16 


0.153 


0.0695 


0.7 


3.95 


2.05 


1.9 


1.08 


0.140 


0.0683 


1.0 


3.85 


1.9 


1.95 


0.97 


0.127 


0.0668 



Table 6: Ring formation function of P,„i, 



^max 
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S/C 


R 


R/S 


2.00 


3.8 


1.5 


2.3 


0.65 


0.090 


0.0600 


2.33 


3.9 


1.8 


2.1 


0.86 


0.120 


0.0667 


2.67 


4.1 


2.2 


1.9 


1.16 


0.153 


0.0695 


3.00 


4.25 


2.5 


1.75 


1.43 


0.190 


0.0760 


3.50 


4.5 


3.2 


1.3 


2.46 


0.253 


0.0791 



Table 7: Ring formation function of a, 



parameter P„,i„ differs from d^^^ in that it is a measure of weighted population density rather 
than of the size of an individual. It is unclear what biological or chemical factors set the 
minimum swarmer biomass, or even if this threshold differs substantially as experimental 
parameters vary. 

As P^in increases, the total time spent swarming, S, decreases, and the time spent in 
consolidation, C increases. This is the same as the response of the system to changes in a,„i„, 
except that the total cycle time, T, remains strongly invariant under changes in d^^^, but 
decreases somewhat as P^i^ increases. The terrace width, R, decreases as P^i^ increases, as 
it does with increases in d^^^. 

The parameter P,„i„ is needed in the models to generate a distinct consolidation phase. 
The loss of a clear distinction between swarming and consolidation as Pj^i„ goes to zero is 
further illustrated in Figure El 
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Figure 3: Swarm-colony front dynamics as a function of Pmin- Shown are curves for the cases 
of Pniin = 0.0, Pniin = 0.2, P^^^ = 0.5, aucl P^i„ = 1.0. The other parameters for these swarm 
colonies are a^^^ = 2.67, Vc = 8, Dq — 2 x 10~^, a^^^ = 0.0, and = 0.5. The curves 
have been realigned in time for better comparison of individual ring formation. This figure 
illustrates the loss of a distinct transition from swarming to consolidation phase as P^,^ goes 
to zero. 
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Table [7| exhibits ring formation as a function of cimax, the dimensionless age at which a 
multinuclear swarmer cell breaks down into mononuclear dividing cells. A swarm colony 
with larger swarmer cells, i.e. larger a^^^^ can be expected to have longer swarming phases, 
5", with wider terrace widths, R. The total cycle time, T, is also seen to increase as a^^^ 
increases, whereas the time spent in consolidation, C, decreases. 

7 Conclusions 

In this paper we have presented models and computational results to understand the periodic 
front dynamics in Proteus mirabilis swarm-colony development. These models differ from 
previous work by placing the most emphasis on known biological mechanisms, in particular 
the local kinetics, or age-structure, and make fewer assumptions about the diffusivity. The 
numerical solutions to the model show the periodicity seen in the biological system, as well 
as the invariance of the cycle time to parameters sensitive to changes in agar and glucose 
concentration, and retain a mechanism for controlling the ratio of swarm to consolidation 
time within an invariant total cycle time. These results were obtained by robust and accurate 
numerical methods designed specifically for age-structured diffusion problems, with proven 
convergence properties. 

Future work on further understanding Proteus swarm-colony development includes ex- 
perimental studies on the forms of the differentiation function, ^. Experiments to determine 
how differentiation from dividing cells to swarmer cells depends on local population density 
would be particularly valuable. 

Other areas of future research are the specifics of how the much larger and broader issue 
of biological motion relates to Proteus motion on an agar surface, and modeling the effects 
of a non-sharp age of septation on the regularity of Proteus front dynamics. Extending the 
models to explicit two-dimensional space to understand the interaction of multiple colonies 
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is also an area of future research. 

Proteus mirabilis should be seen as a model system with relevance to other biological 
problems. The further mathematical understanding of age and space structured models - 
and the numerics to solve them - constitutes an important foundational area of research. 
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A Numerical Solutions of the Esipov-Shapiro Model 

The Esipov-Shapiro model, as described by the relevant equations and parameters in [T6] . 
has no spatial dynamics whatsoever. 

We make several natural changes in the Esipov-Shapiro model to generate something 
resembling Proteus mirabilis swarm-colony development. These changes, although not con- 
tained in their equations, are suggested in the text. 

We set the units of the parameter Dq to cm^/hr rather than cm^/sec. Use of the latter 
causes very rapid swarming that covers the entire domain in a short time. It is reasonable 
to believe that this is a typo, as all other parameters in use hours as the time unit. The 
second change is the use of swarmer biomass rather than numerical density as the relevant 
measure of total population. Although the formal mathematical definition of P, (page 255 
of [inj) is given by Pg = dO Ps(r, 6, t), the statement at the top of page 253, 

Surface densities of swimmer and swarmer cell populations will be indicated 
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by the capital letters Pc and Ps, respectively (c for consolidation phase, s for 
swarming phase). While we do not resolve swimmers in age, the surface density 
of swarmers is related to a biomass weighted average over the ages, with e^/"^"* 
being the contribution of age 6. 

along with personal communication with Sergei Esipov (July, 1997), and the statement by 
the authors in j44j that their models in |16| use an "age- weighted swarmer cell density", 
indicates that a mass-weighted Pg was used in their numerical computations. This is a 
necessary reinterpretation. The use of a non-weighted Ps precludes the upper trigger of the 
memory field needed to instigate swarming, Ps,max, from ever being reached, which in turn 
results in no swarm-colony development whatsoever. 

Another change in the Esipov-Shapiro model needed to get some sort of swarming and 
consolidation behavior is the incorporation of what they refer to as the "drivers license" age, 
^miii' explicitly into the diffusivity. As the equations are written, 9jj^[^ enters the system 
implicitly only through the definition of P^. Thus, small swarmers move with larger swarmers 
in the "rafts" ; they just do not contribute to a raft's motility. This is a reasonable assumption 
- nothing prevents them from being caught up in the flow - and one made in the models 
presented in this paper. But because we use a lag phase in septation, the presence or absence 
of "underage" swarmers in rafts does not qualitatively affect our model. However, since it 
does not contain such a lag, the Esipov-Shapiro model does not develop a consohdation 
phase after the first swarm phase, but rather the swarmers develop a self-sustaining soliton 
caused by swarmers of maximum age that dedifferentiate into swimmers (dividing cells,) 
immediately differentiating into new swarmers. Although this does not reflect the observed 
behavior of Proteus, it is an interesting mathematical object, which we display in Figure |3] 

If we reinterpret the meaning of the "drivers license" age to be that small swarmers 
do not move at all, we need to multiply the diffusivity, D, as written in Equation (4a) on 
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Figure 4: Profiles of swarmer biomass for the Esipov-Sliapiro model without OiOmin) 
diffusivity. The model does not develop a consolidation phase after the first swarm phase, 
but rather the swarmers develop a self-sustaining soliton where swarmers of maximum age 
that dedifferentiate into dividing cells immediately differentiate into new swarmers. The 
parameters used in this computation are taken from the caption of Figure 4 in jTHj, with the 
choice of = 0.045. 

page 256 of jTHI, by the Heaviside function QiOiain)- Doing so results in a system with a 
swarming/ consolidation cycle with fixed time lengths, but not the regularly spaced terraces 
seen in the experiments. The terrace widths decrease in length each cycle. This is illustrated 
for several values of Dq in the Esipov- Shapiro model in Figure |S1 

The cause of the spatial regularity seen in the computations in jTH] is quite possibly due 
to computational effects. The numerical integration in age and time in [T^ is not done 
along characteristics, but rather consists of transforming age to a stage-structured variable 
dependent upon discretization parameters, adding regularity to a system that is not really 
there. Another possible source of error is that the space and time integration, as written in 
Equations (6b)- (6d) on page 257 of has the property that the colony front can move at 
most one spatial interval per time step. Thus, inaccurate numerics can change the nature of 
the bacterial front. 
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Figure 5: Swarm-colony front dynamics of the Esipov-Shapiro model as a function of their 
Dq. The parameters used in these computations are taken from the caption of Figure 4 in 
nil, including the choice of Dq = (0.045, 0.036, 0.015). 



We chose a smooth piecewise cubic function for the form of ^ in Equation ()4.1i|) for numerical 
reasons. In this appendix, we argue that it is the lag in ^, not the shape of the curve, that 
is relevant in modeling Proteus swarm-colony development. We present two sets of results 
from using two different piecewise constant functions with the same area under the curve as 
the piecewise cubic. The first, "fat" one, 



has the same support, but shorter maximum height, than the cubic. The second, "skinny" 
one. 



B Use of Different Forms of ^ 




(B.5) 
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2.1 
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0.134 



Table 8: Ring formation as a function of Vc using Equation ()B.5|) . 
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3.2 


2.2 
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0.117 
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3.5 


2.3 


1.2 


0.147 
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4.0 


2.3 


1.7 


0.164 


16 


4.5 


2.3 


2.2 


0.176 



Table 9: Ring formation as a function of Vc using Equation ()B.6|) . 

has the same maximum height, but narrower support, than the cubic. Using Equation ()B.5jl 
instead of Equation ()4.1i|) gives Table |H1 in place of the same data in Tabled Using Equation 
(1B.6|1 instead of Equation ()4.1i|l gives Table El in place of those data in Table IH 

The temporal and spatial regularity seen in Proteus and the computations presented 
earlier with a piecewise cubic ^ are also present in these computations. These results indicate 
that changing the form of ^ does not alter the qualitative behavior of the system and has 
only a minor quantitative effect on the solutions. Nonetheless, the shape of ^, and not only 
its mass and support, is an interesting topic for experiments. 
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